#----------------------
#Likelihood
#----------------------
likelihood <- function (x, dataX, dataY1, U = NULL, center = NULL, 
	pred = FALSE, uniforming = 0, raw_val = TRUE)
{
	if (is.null (U)) {
		xt <- rep (0, length (x))
		#PARAMETERS
		xt <- x
	} else {
		xt <- Solve (U, x)
		xt <- xt + center  
	}
	y1_pred <- log.model (xt, dataX)
	sigma1 <- exp (xt [length (x)])#precision error is specified?
	error1 <- (log (dataY1) - y1_pred)
	prob_Y1 <- dnorm (as.matrix (error1), mean = 0, sd = sigma1, 
		log = TRUE)
	if (raw_val == TRUE) {
		if (pred == FALSE) {
			c (xt, (uniforming + sum (prob_Y1)))
		} else {
      		list(sum (prob_Y1), y1_pred, xt, sigma1)
		}               
	} else {
		if (pred == FALSE) {
			c (x, (uniforming + sum (prob_Y1)))
		} else {
      		list (sum (prob_Y1), y1_pred, x, sigma1)
		}
	}
}
#-------------------------
#Prior distribution
#-------------------------
prior <- function (x, m, variance)
{
	if (x[3] > x[1]) {
#     	dmvnorm(x, m, variance, log = TRUE)
      	0
	} else {
		0
	}
}
#----------------------
#Posterior
#----------------------
joint <- function (x, dataX, dataY1, U = NULL, center = NULL, uniforming = 0,
	raw_val = TRUE, logf = TRUE,  ...)
{
	if (logf) {
      	prior (x, ...) + likelihood (x, dataX, dataY1, U, center, 
			pred = FALSE, uniforming, raw_val)
	} else {
		prior (x, ...) * likelihood (x, dataX, dataY1, U, center, 
			pred = FALSE, uniforming, raw_val)
	}
}

#---------------------------------
#Obtaining q-points from range
#---------------------------------
legendre_conversion <- function (var_range, q_points)
{
	upper <- var_range [2]
	lower <- var_range [1]
	eval_points_var <- (upper + lower) / 2 + (upper - lower) / 2 * q_points
	eval_points_var
} 

#-----------------
#Range definition
#-----------------
range_definition <- function (eval_var, var_range, means, q_points, dataX, 
	dataY, threshold, U = NULL, center = NULL, raw_val = TRUE)
{
	#defining evaluating points from quad. points
	eval_points_var <- legendre_conversion (var_range, q_points)
	#creating a matrix with mean values
	eval_points <- matrix (data <- rep (means, length (q_points)), 
		ncol = length(means), nrow = length (q_points), byrow = TRUE) 
	#subst. values of the var to evaluate
	eval_points [, eval_var] <- eval_points_var 
	results <- apply (eval_points, MARGIN = 1, 
		function (x)
		{
			joint (x, dataX, dataY, U, center, uniforming = 0,
				raw_val, logf = TRUE)
		}) [(length (means) + 1), ]
	significant <- (results >= threshold) #comparing with threshold 
	#"good" points position
	positions <- seq (length (eval_points_var)) [significant] 
	min_value <- min	(positions) #min "good" value position
	max_value <- max	(positions) #max "good" value position
	min_a <- if (min_value != 1)
		{
			min_value - 1
		} else {
			min_value
		}# is it the given limit?
	max_b <- if (max_value != length (q_points))
		{
			max_value + 1
		} else {
			max_value
		} #given lim.?
	matrix (data <- c (eval_points_var [min_a], eval_points_var [min_value], 
		eval_points_var [max_value], eval_points_var [max_b]), ncol = 2, 
		nrow = 2, byrow = T) 
}                